3. Motifs

In this tutorial, waveform motif identification will be demonstrated with simulated and ECoG data.

Simulated Data

First, a timeseries containing two bursty, asymmetrical sine waves with a 1/f background will be simulated.

import numpy as np

from neurodsp.sim import sim_variable_oscillation, sim_powerlaw
from neurodsp.plts import plot_time_series, plot_power_spectra
from neurodsp.spectral import compute_spectrum
from neurodsp.filt import filter_signal
from neurodsp.utils.norm import normalize_sig

from fooof import FOOOF

from ndspflow.motif import Motif

# Settings
fs = 1000
freq = 20
n_cycles = 5
n_choices = 80
n_seconds = 1/freq * n_cycles * n_choices

# Define signals
sig_lo_rdsym = sim_variable_oscillation(None, fs, freq, cycle='asine',
                                        rdsym=[.2] * n_cycles, phase='min')
sig_hi_rdsym = sim_variable_oscillation(None, fs, freq, cycle='asine',
                                        rdsym=[.8] * n_cycles, phase='min')
sig_ap = sim_powerlaw(n_seconds, fs, exponent=-2)
sig = np.zeros_like(sig_ap)

# Sample defintions
cyc_len = len(sig_lo_rdsym)
cyc_starts = np.arange(0, len(sig), cyc_len)
cyc_ends = np.arange(cyc_len, len(sig)+cyc_len, cyc_len)

# Randomly choose between .2 rdsym, .8 rdsym, or 1/f powerlaw
choices = [np.zeros_like(sig_lo_rdsym), sig_lo_rdsym, sig_hi_rdsym]
for start, end in zip(cyc_starts, cyc_ends):
    choice = np.random.choice([0, 1, 2])
    sig[start:end] = np.sum((choices[choice] * 0.5,
                             sig_ap[start:end]), axis=0)

sig_pe = sig - sig_ap

# Plot simulated signal
times = np.arange(0, len(sig)/fs, 1/fs)

plot_time_series(times, sig, title='Combined Signal', xlim=(0, 20))
plot_time_series(times, sig_pe, title='Periodic Component', xlim=(0, 20))
plot_time_series(times, sig_ap, title='Aperiodic Component', xlim=(0, 20))

Spectral Parametrization

Identifying motifs requires spectral parameterization to define center frequencies and bandwidths. Note, a list of tuples, defining center frequencies and bandwidth may be alternatively used.

# Spectral parametrization
f_range = (1, 100)
freqs, powers = compute_spectrum(sig, fs, f_range=f_range)

fm = FOOOF(peak_width_limits=(2, 10), peak_threshold=2.5)
fm.fit(freqs, powers, freq_range=f_range)

Motif Extraction

Motif extraction involves using Bycyles to segement cycles by defining extrema and zero-crossings. Cycles within a peaks bandwidth are resampled to the center freqeuency, allowing a mean motif to be found. Each segmented waveform found within the bandwidth of a spectral peak is resample to the center frequency, allowing a mean waveform to be associated with each peak. Multiple waveforms at a given peak are distinguished using k-means clustering. Cycles highly correlated (r >= 0.8) with a motif are extracted and a refined motif is recomputed.

# Motif extraction
motif = Motif()
motif.fit(fm, sig, fs)
motif.plot()

Signal Decomposition

Once motifs are found, they are affine transformed to each individual cycles. The resulting transformation matrix describes cycle-by-cycle variation. Transformed motifs provide an estimate of the periodic signal. The difference between transformed motifs and cycle waveforms provide an estimate of the aperiodic signal. Spectra of the decomposed signals may be used to validate the aperiodic/periodic separation in the frequency domain.

# Signal decomposition
motif.decompose()
motif.plot_decompose(0)

# Plot spectra of the decomposed signals
motif.plot_spectra(0)

ECoG Data

Motif extraction and signal decomposition will be applied to a pubically available EcoG recording.

Fedele, T., Boran, E., Chirkov, V., Hilfiker, P., Grunwald, T., Stieglitz, L., … Sarnthein, J. (2021). Dataset of spiking and LFP activity invasively recorded in the human amygdala during aversive dynamic stimuli. Scientific Data, 8(1). doi:10.1038/s41597-020-00790-x

# Normalize and filter the data
f_range = (1, 100)
sig_ecog = np.load('data/ecog.npy')
fs = 2000

sig_ecog = normalize_sig(sig_ecog, mean=0, variance=1)
sig_ecog = filter_signal(sig_ecog, fs, 'bandpass', f_range, remove_edges=False)

# Spectral parameterization
freqs, powers = compute_spectrum(sig_ecog, fs, f_range=f_range)
fm_ecog = FOOOF(peak_width_limits=(2, 8), peak_threshold=2, max_n_peaks=3, aperiodic_mode='knee')
fm_ecog.fit(freqs, powers, f_range)

# Motif extraction
motif_ecog = Motif(min_clust_score=.2, corr_thresh=.1, max_clusters=4)
motif_ecog.fit(fm_ecog, sig_ecog, fs)
motif_ecog.plot(plot_fm_kwargs={'log_freqs':True})

# Signal Decomposition
motif_ecog.decompose()
motif_ecog.plot_decompose(0, xlim=(20, 24))

# Plot decomposed spectra
motif_ecog.plot_spectra(0)

# Plot the affine transformation parameters
motif_ecog.plot_transform(0, xlim=(20, 24))